Section: New Results

Incremental algorithms for long-range interactions

Participants : Semeho Edorh, Stephane Redon.

Numerical simulations of molecular dynamics (MD) are very expensive in terms of CPU resources. During Molecular dynamics simulations, the most CPU intensive task is the evaluation of the interaction potential [78] . Due to the large number of particles involved, updating this potential may have, at each time-step, a very high computational cost.

In large crystalline ionic system, Ewald summation is the most popular method for computing electrostatic interactions. It rewrites the interaction potential $\phi$ as the sum of a short-range term and a long-range term. Ewald summation using optimal parameters requires $𝒪\left(N^{3/2}\right)$ operations [47] , [30] but it can be modified so that it involves only about $𝒪(NlogN)$ operations [31] , [85] by using the Fast Fourier Transform.

We want to develop a new approach that can reduce the computational cost by using incremental algorithms. The key-idea is to use, for each time-step of the simulation, information that we have computed in previous steps.

The Particle Mesh Ewald (PME) algorithm developped by Darden et al. is the most successful approach for computing long range interactions. In the particle mesh method, just as in standard Ewald summation, the generic interaction potential is separated into two terms.The so-called short-range contribution can be easily calculated in a direct space by using truncation methods. Where as the long-range contribution is calculated using two Fast Fourier transforms ($Nlog\left(N\right)$ algorithm ). In pratice, the long range contribution algorithm boil than to [30] :

• Map particle charge density $Q$ to a mesh

• Compute the forward Fast Fourier Transform of the approximation $Qm$ of charge density on the mesh

• Multiply $Qm$ by a green function (related to the choice of the mesh).

• Compute a backward Fast Fourier Transform of the result.

• Retrieve the long-range contribution potential by interpolating the previous result at particles positions.

We modified this algorithm to make it incremental. We started from the PME implementation in LAMMPS. Instead of mapping the charge density to the mesh, we mapped the increment of density $dQ$ to the mesh. The FFT solver KissFFT is based on a divide-and-conquer algorithm. We built a sparse input solver as a modified version of FFT solver which computes only needed (non trivial) operations [74] . We built also a sparse output solver inspired by the algorithm proposed by Katabi et al. [36] . Unfortunately, we did not get significant speed-ups with these modifications.

We decided to compute the increment of the long-range contribution related to the increment of density $dQ$ by using multi-resolution methods. These methods are slower than PME but have better adaptive behavior. The multigrid approach was chosen because of its $𝒪\left(N\right)$ behavior and its good scalability [13] . We are currently developing an adaptive multigrid method.